logo

1 Load R package

Currently only on github - but will be available on CRAN soon

# remove.packages("brainMapR")
devtools::install_github("baptisteCD/brainMapR")
library(brainMapR)

Some brain association maps (for age at MRI) are given with the package. They are stored in the package source folder and you can find the paths using the commands below.

system.file("extdata/07_BWAS_MOA_realPheno", "", package = "brainMapR",mustWork = TRUE)

system.file("extdata/03_BWAS_uncorrected_realPheno", "", package = "brainMapR",mustWork = TRUE)

system.file("extdata/03_BWAS_uncorrected_realPheno", "BWAS_Age_MRI.linear", package = "brainMapR",mustWork = TRUE)

2 Extract cluster characteristics

for (pheno in c("Age_MRI","body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "sexuses_datacoding_9_f31_0_0","smoking_status_f20116_2_0", "sexuses_datacoding_9_f31_0_0") ){
for (scenario in c("03_BWAS_uncorrected_realPheno", "09_BWAS_ICVagesex_realPheno", "04_BWAS_5globalPCs_realPheno", "05_BWAS_10globalPCs_realPheno", "06_BWAS_10specificPCs_realPheno", "07_BWAS_MOA_realPheno",  "10_BWAS_MOA_multiORM_QC_realPheno", "15_BWAS_MOA_allModa_FE_realPheno")){
  
identifyClustersBWAS(pathToTheBWASresults =system.file(paste0("extdata/", scenario), "", package = "brainMapR",mustWork = TRUE), bwasFile = paste0("BWAS_",pheno), outputFolder = "../", signifThreshold = 5e-8)
  print(paste(pheno, scenario))
}
}

3 Create Legend bar

Legend to be used in plots

cols=c(RColorBrewer::brewer.pal(n = 10, name = "RdYlBu")[10:6],RColorBrewer::brewer.pal(n = 10, name = "RdYlBu")[5:1]) # Select palette colours

png("legendbar.png", width = 5, height = 15, units = "cm", res = 400)
par(mar=c(2,4,2,2))
my.colors = colorRampPalette(cols)
z=matrix(1:100,nrow=1)
x=1
y=seq(-0.4,0.4,len=100) # supposing 3 and 2345 are the range of your data
image(x,y,z,col=my.colors(20),axes=FALSE,xlab="",ylab="", main="Correlation")
axis(2)
dev.off()
library(knitr)
include_graphics(path = "examplePlots/legendbar.png", dpi = 50)

4 Make brain plots

Brain surface plots, made in R using the rgl package. The cortical surface can take a while to render but the plotting is very versatile.
The brain plot functions use the summary files created using identifyClustersBWAS().
Brain plots will silently open the phenotype file (located in UKB_phenotypes15K) to scale association betas into correlation coefficients.

library(viridis)
cols=viridis_pal(option = "C")(7)

for (var in c("Age_MRI","body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "sexuses_datacoding_9_f31_0_0","smoking_status_f20116_2_0", "sexuses_datacoding_9_f31_0_0") ){
  for (scenario in c("03_BWAS_uncorrected_realPheno", "09_BWAS_ICVagesex_realPheno", "04_BWAS_5globalPCs_realPheno", "05_BWAS_10globalPCs_realPheno", "06_BWAS_10specificPCs_realPheno", "07_BWAS_MOA_realPheno",  "10_BWAS_MOA_multiORM_QC_realPheno", "15_BWAS_MOA_allModa_FE_realPheno")){

plotCortical(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")

plotSubcortical(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")

plotSubcortical_flat(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")

combineCorticalSubcorticalPlots(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../" )
}
}

4.1 Cortical panels

Left cortical surface area associations with age - GLM without covariates

library(knitr)
hemi="lh"
moda="area"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))

> Right cortical surface area associations with age - GLM without covariates

library(knitr)
hemi="rh"
moda="area"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))

Left cortical thickness associations with age - GLM without covariates

library(knitr)
hemi="lh"
moda="thickness"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))

Right cortical thickness associations with age - GLM without covariates

library(knitr)
hemi="rh"
moda="thickness"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))

4.2 Subcortical plots

Sucbcortical surface area associations with age - GLM without covariates

library(knitr)
hemi="rh"
moda="LogJacs"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))

hemi="lh"
moda="LogJacs"
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_inside.png"))
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/BWAS_Age_MRI_", hemi, "_", moda, "_clustersAndCoordinates_outside.png"))

4.3 Combined plot

include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/Plots_CombinedAge_MRI.png"))

5 Produce Manhattan plots

library(viridis)
cols=viridis_pal(option = "C")(7)

BrainMapManhattanPlot(inputPath = system.file("extdata/07_BWAS_MOA_realPheno", "", package = "brainMapR",mustWork = TRUE), path = "03_BWAS_uncorrected_realPheno/", bwasFile = "BWAS_Age_MRI.linear", yMax = 120, signifThreshold =5e-8,  phenotypeLabel = "Age")
include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/Manhathan_BWAS_Age_MRI.linear_simple.png"))

6 Table of associated regions and description

UKB_BWASvars_labels.txt contains 3 columns: variable names, category (e.g. cognition, lifestyle) and variable label (used for plots and tables)

pheno=read.table("UKB_BWASvars_labels.txt", header=T, stringsAsFactors = F)
pheno=pheno[which(pheno$variable %in% c( "body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "Age_MRI","sexuses_datacoding_9_f31_0_0", "smoking_status_f20116_2_0")),]

vall=NULL
for (iii in 1:5){
  all=NULL
  for (scenario in c("03_BWAS_uncorrected_realPheno","09_BWAS_ICVagesex_realPheno", "05_BWAS_10globalPCs_realPheno", "07_BWAS_MOA_realPheno", "15_BWAS_MOA_allModa_FE_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno")){
  
  vres=read.csv(paste0("../../",scenario, "/", "BWAS_signif_", pheno$variable_clean[iii], ".csv" ))
  res=read.table(paste0("../../",scenario, "/", "Results_clusterFWER_", pheno$variable[iii], ".txt" ), header=T)
  nclust=res[,c("NumberClusters_thickness", "NumberClusters_area", "NumberClusters_thick", "NumberClusters_LogJacs")]
  maxclust=res[,c("maxSizeCluster_thickness", "maxSizeCluster_area", "maxSizeCluster_thick", "maxSizeCluster_LogJacs")]
  
  all=rbind(all,  dim(vres)[1], sum(nclust), max(maxclust, na.rm = T))

  }
  vall=cbind(vall, all)
}

colnames(vall)=pheno$variable_clean
writexl::write_xlsx(as.data.frame(vall), "../../Summary_BWAS_realPheno_bonferroniNtraits_R1IEEE.xlsx")


# N clusters cortex
vall=NULL
for (iii in c(1,2,3,5)){
  all=NULL
  for (scenario in c("03_BWAS_uncorrected_realPheno","09_BWAS_ICVagesex_realPheno", "05_BWAS_10globalPCs_realPheno", "07_BWAS_MOA_realPheno","15_BWAS_MOA_allModa_FE_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno")){
  #vres=read.csv(paste0("../../",scenario, "/", "BWAS_signif_", pheno$variable_clean[iii], ".csv" ))
  res=read.table(paste0("../../",scenario, "/", "Results_clusterFWER_", pheno$variable[iii], ".txt" ), header=T)
  nclust=res[,c("NumberClusters_thickness", "NumberClusters_area", "NumberClusters_thick", "NumberClusters_LogJacs")]
  maxclust=res[,c("maxSizeCluster_thickness", "maxSizeCluster_area", "maxSizeCluster_thick", "maxSizeCluster_LogJacs")]
  all=rbind(all, sum(nclust[,1:2]))
  }
  vall=cbind(vall, all)
}
colnames(vall)=pheno$variable_clean



 



A work by by [Baptiste Couvy-Duchesne] - 08 October 2021